NASA Technical Memorandum 101002 US AAVSCOM Technical Report 88-A-006 


NAS A-TM- 1 0 1 002 19880017325 


A General-Purpose Approach to 
Computer-Aided Dynamic 
Analysis of a Flexible Helicopter 

Om P. Agrawal 


For Reference 


nor to ts IAJ 2 H rcou rac r.oou 


July 1988 


rj / * r. F r. 


NASA 

National Aeronautics and 
Space Administration 



SYSTEMS COMMAND 


AVIATION RESEARCH AND 
TECHNOLOGY ACTIVITY 




National Aeronautics and 
Space Administration 

Ames Research Center 

Moffett Field, California 94035 

September 16, 1988 


Reply to Attn of: Shelley Scarich 


ERRATUM NASA TM 101002 

Please replace the copy (or copies) you were sent of the attached TM by 
Om P. Agrawal with these copies. Copies were mailed out about 9/14/88. 
If you had not received a copy, this/these is/are the file copy(ies) . 
Thanks . 





NASA 



NASA Technical Memorandum 101002 


USAAVSCOM Technical Report 88-A-006 


A General-Purpose Approach to 
Computer-Aided Dynamic 
Analysis of a Flexible Helicopter 

Om P. Agrawal, Southern Illinois University at Carbondale, Carbondale, lllionois 


July 1988 


NASA 

National Aeronautics and 
Space Administration 



SYSTEMS COMMAND 


Ames Research Center 

Moffett Reid, California 94035 


AVIATION RESEARCH AND 
TECHNOLOGY ACTIVITY 
MOFFETT FIELD, CA 94305-1099 







SUMMARY 


This report describes a general-purpose mathematical model for dynamic analysis of a helicopter con- 
sisting of flexible and/or rigid bodies (components) that may undergo large translations and large rotations. 
In this formulation, two sets of generalized coordinates are used, namely rigid-body and elastic. The rigid 
body coordinates define the location and the orientation of a body coordinate frame (global frame) with 
respect to an inertial frame. The elastic coordinates are introduced using a finite element approach in order 
to model flexible components with complex shapes. The compatibility conditions between two adjacent 
elements of a flexible body are imposed using a Boolean matrix, whereas the compatibility conditions be- 
tween two adjacent bodies are imposed using the Lagrange multiplier approach. The virtual work approach 
is used to derive the equations of motion. This leads to a system of differential and algebraic equations 
that are highly nonlinear due to large rotations of the components of the system. Some numerical schemes 
to solve the resulting set of equations are also discussed. 

The formulation presented is general. When compared with a hard-wired formulation, this approach 
has several advantages. 1) The flexible model of each component in the system is processed identically. 
Therefore, only one set of subroutines is required to generate several matrices such as inertia and stiffness 
matrices. 2) The constant submatrices associated with a component remain constant throughout the sim- 
ulation and the sparsity of these component matrices is preserved because the matrices associated with a 
component do not depend on the generalized coordinates of the other components. 3) The form of con- 
straint equations depends upon the type of kinematic joint and involves the generalized coordinates of 
the two components participating in the kinematic joint. Therefore, a library of constraint elements can 
be developed and the kinematic constraint between two adjacent bodies can be imposed in an automated 
fashion. 4) This formulation is capable of generating equations of motion for a system with any number 
of flexible and/or rigid bodies connected in an arbitrary fashion. 5) The Lagrange multipliers which are a 
byproduct of this formulation provide the reaction forces and torques at the joints. 
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1. INTRODUCTION 


Analysis of a helicopter consisting of interconnected rigid and/or flexible components, each of which 
may undergo large translations and large rotations, requires a mathematical model that accounts for nonlin- 
ear kinematic constraints, nonlinear couplings between rigid motion and elastic deformations, and inertia 
variations resulting from elastic deformations. In mathematical models developed in the past (refs. 1-6), 
the nonlinear kinematic constraints were removed explicitly and the equations of motion were developed 
in terms of a minimum number of generalized coordinates. Such models are, although quite appealing 
because they lead to a set of minimum number of differential equations, quite involved in respect of their 
mathematical formulations. In contrast to automated coupling, these models are often limited in their 
flexibility to represent every possible kinematic joint in terms of the minimum number of generalized co- 
ordinates considered, and are also limited in their expandability to a new system since they are based on a 
specific topology of the system. Furthermore, these models destroy the sparsity and the constant matrices 
appearing in the component equations of motion, which makes it difficult to take advantage of the sparsity 
and modal analysis for space economy and dimension reduction. Finally, these models often introduce 
singularities in the models which require special care. 

It is clear that a new mathematical model for analysis of a helicopter is desired. This report develops 
a general purpose model for analysis of a helicopter and identifies features of this model which makes it 
superior to the previously developed models. 

2. GENERALIZED COORDINATES AND VIRTUAL WORK DUE TO INERTIA FORCES 


In order to specify the configuration of a body (component) of a helicopter system, it is necessary to 
define a set of generalized coordinates such that the absolute position of a point in the body can be defined in 
terms of those generalized coordinates. As shown in figure 1, XYZ is the inertial frame and X'Y'Z' is the 
coordinate frame (or the global frame) associated with body i. The body i may be a fuselage, a rotor, or any 
other component of the helicopter system. In the discussion to follow, the coordinate system associated 
with body i will be called frame i. Using figure 1, the position vector of point p,r p , can be written as 
(refs. 7, 8), 


r p = R i + A i s (2.1) 

where R'( = [a*, y\ 2 *]) is the position vector of the origin of frame i with respect to the inertial frame, A' 
is the orthonormal transformation matrix from frame i to the inertial frame, and s is the position vector of 
point p with respect to frame i. Vector s may be written as 


s = so + u (2 .2) 

where s 0 is the undeformed location of point p with respect to frame i, and u is the vector field of elastic 
deformations. Body i is idealized as a finite collection of structural elements having nn number of nodes 
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and rirf number of nodal degrees of freedom. Assuming elastic deformations to be small, the vector field u 
can be approximated as 


u = N i q i (2.3) 

where q x is a m x 1 vector of nodal elastic deformations of body i and N ‘ is the shape function. Note that 
this model allows one to have different numbers of degrees of freedom at each node. Although the virtual 
work expression due to inertia forces (and subsequently the equations of motion) can be developed in terms 
of field variable u, the finite element approach is used here to allow modeling of flexible components with 
complex shapes. 

Substituting equations (2.2) and (2.3) into equation (2.1), the position vector of point p can be approx- 
imated as 


(2.4a) 
(2.46) 

Differentiating equation (2.4) with respect to time, the velocity vector of point p is obtained. 


r p = R' + A'(s 0 + u) 

= R' + A'( s 0 + JV V) 


r p =R i + A i s + A i i i (2.5 a) 

= R' + A's + A'N'q' (2.56) 

where (•) denotes differentiation with respect to time. The second term in the right-hand side of equation 
(2.5) can be written as 


1 8A * 

A i s = ^—.9 } s = B i (9 i ,s)e< 


U d6 > 


( 2 . 6 ) 


where 6j(j = 1 , ..., 1) are the rotational degrees of freedom of the ith global frame, and 6' is the vector of 
e jU = l ( ...,/),i.e., 


e i 


0 ‘ = 


0 / 


(2.7) 


These rotational degrees of freedom can be Euler angles, Euler-like angles, Euler parameters, or other 
rotational coordinates including quasi-rotational coordinates. In most cases 1=3. However, in the case 
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of Euler parameters, 1=4. Selection of rotational degrees of freedom has been an issue of much debate. 
For small rotations, most rotational degrees of freedom do not present any problem. For large rotations, 
however, certain rotational degrees of freedom may lead to singularity problems. 

Since the matrix B * is linear in its second argument, it follows that 


B\9,u + v) = B\6,u) + B\6,v) 


(2.8) 


This property of matrix B x is used later to isolate the rigid body, the elastic body, and the coupling terms. 
The matrix B' is linear in its first argument if Euler parameters are used to define the rotational degrees of 
freedom. Using equation (2.6), the velocity vector of point p can be expressed as 


f p = R* + B\ 6, s) 9 { + A'ii (2.9 a) 

= ti + B\9, s) 0‘‘ + A'N*? (2.9 6) 


Define 


P< = 


' R { ‘ 
6 i 


( 2 . 10 ) 


to be the vector of generalized coordinates of body i. The first two components of P l represent rigid body 
degrees of freedom and q' represents the elastic coordinates. Using equation (2.10), the velocity of point 
p becomes 


f p =[l B\6,a) A i N i ]P i (2.11) 

where I is a 3x3 identity matrix. Differentiating equation (2.11) with respect to time, one obtains the 
acceleration vector of point p as 

r p =[ I BKO, s ) A i N i ] P i + B\ 6 , s) 0*' + AWq* (2.12 a) 

= [/ B { (9,s) A*N' ] P* + B\ 9, s) 6' + B\ 9, it) 0* (2.126) 


Similarly the virtual displacement of point p can be given by 
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8r p = SR* + 8A*s + A { 8u 


(2.13o) 


= 8R' + B'(6,s)86 i + A'N'Sq' 


(2.13 b) 


where 8( ) represents the variation of a parameter. 

The virtual work due to inertia forces can now be written as 



(2.146) 


where dm * is the differential mass of a small volume at point p, m* is the total mass of body i, and the 
superscript T represents the transpose of a vector (or a matrix). Note that the virtual work due to inertia 
forces contains three distinct virtual work expressions, namely, virtual work due to: (1) virtual translations 
of frame i, (2) virtual rotations of frame i, and (3) virtual elastic deformations. Substituting equations 
(2.12) and (2.13) into equation (2.14) and using equation (2.10), the three virtual work expressions due to 
inertia forces can be written as 

1) Virtual work due to virtual translation of frame i: 

= SB? [ [ I B'(O.s) 

J TO* 

+ 8R? [ .(£‘(0, s)6 i + B‘(0, u)^) dm*' '(2.15a) 

Jm i 


2) Virtual work due to virtual rotation of frame i: 


8W'r 


I(rot) 


= 86? [ [ B?(6, s) B iT (6, s)B i (9, s ) £‘ T (0, s)A i N i ] P* dm * 
J m i *■ J 

+ 86? f (B?(6, s)B i (6,s)6 i + B?(6,s)B i (6,u)6 i ) dm *' 


(2.156) 


3) Virtual work due to virtual elastic deformation: 


8Wi (elas) = 8q? J . [ N iT A? N? A? B i (6, s ) ] P’ dm { 

+ 8q? [ ( N iT A iT B i (6 > s ) 6* + N? A? B\6, u) P) dm { (2.15c) 

Jm { 
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Combining equations (2.15a) to (2.15c), the expression for virtual work due to inertia forces of body i 
becomes 


6W I = 6P f M i P i + 8P {r F i (2.16) 

where M* is the mass matrix of body i and F' is a vector of nonlinear terms that includes generalized forces 
due to centrifugal and Coriolis accelerations. Matrix M l and vector F' are, respectively, given by 


Kr 

M' tB 

M' rq 


Mle 

m q 

L< 

•T 

M % 

K 


(2.17) 


and 


where 



BK0,s) + B*'(0,u) 
B' r (9,s)B i (9,s ) + B' T (9,s)B'(9,u) 
N iT A' T B i (9,s) + N^A^B^d,^) 


9' dm' 


(2.18) 



I dm' = m'l 


is the diagonal rigid-body translational mass matrix, 


M' BB = [ B ,r (9, s)B\9,s) dm { 

Jm * 


is the rotational inertia matrix, 


M' = [ N iT N i dm i 
™ Jm * 


is the elastic mass matrix, 


ML= [ B'(9,s)dm' 

Jm' 


(2.19) 


( 2 . 20 ) 


( 2 . 21 ) 


( 2 . 22 ) 
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(2.23) 


M* = f A'N'dm' = A' f N i dm i 
H Jm' Jm i 


and 


ML = [ B' T (6, s) A'N'dm' 

Jm x 


(2.24) 


are inertia coupling matrices. 

The matrices that appear in equation (2.17) (and equations (2.19) to (2.24)) need further examination. 
Matrices M* r and M qq are constant matrices and therefore are generated only once. Matrix M qq (defined in 
eq. (2.21)) is the same as the mass matrices which appear in the usual finite element formulation of struc- 
tural systems. In contrast to the rigid body matrices, matrices M*g and M' g$ depend on elastic deformations 
of the body. Using equations (2.2) and (2.8), matrix M*g can be written as 


M' t6 = f B\9, so) dm' + [ B\6,u) dm* 
Jm 1 Jm 1 


(2.25) 


where the first term on the right hand side is associated with the moment of mass of the undeformed body 
about its body frame axes and the second term is associated with the moment of mass about the body frame 
axes due to elastic deformations. The first term is zero if the origin of the undeformed body frame is taken 
at the center of mass of the body. 

Using equations (2.2) and (2.8) once again, the rotational inertia matrix Mg 0 can be written as 
Mgg = [ B iT (6,s 0 )B i (e,so)dm i + [ (B' T (6, s 0 )B i (e, u) 

Jm 1 Jm % 

+ B ir (6,u)B i (6, so) + B {r (6,u)B i (6,u))dm i (2.26) 

where the first term on the right hand side is the rotational inertia matrix associated with rigid body rotation 
and the other terms describe contributions due to elastic deformations. Using equation (2.6), it is clear that 
the shape function and the rotational degrees of freedom in equations (2.22) and (2.23) are isolated from 
each other and therefore the integral 


/ N'dm* (2.27) 

Jm' 

is generated only once. Similarly, the rotational degrees of freedom and the elastic deformation in equations 
(2.18), (2.20), and (2.24) can be isolated. The mathematical manipulations in this case depend on the 
rotational degrees of freedom selected. 
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From the above discussion, it is clear that the matrix M' is highly nonlinear because of large rotations 
and elastic deformations. Furthermore, rigid body motion and elastic deformations are strongly coupled. 
An accurate model, which allows the components of the helicopter system to undergo large rotations, must 
include these nonlinear and inertia coupling effects. Neglecting these effects in helicopter dynamics may 
result in a significant error in the dynamic response of the helicopter since generalized forces are quite 
sensitive to generalized coordinates. 

The virtual work expression due to inertia forces presented above is general and applies to each body 
equally. Therefore, the total virtual work due to inertia forces is 


NB 

SW I = J 2 5w i (2.28) 

»=1 

where NB is the number of components or bodies in the system. The generalized coordinates of each 
component need not be independent. If one needs to develop the virtual work 8Wi in terms of independent 
generalized coordinates, then equation (2.28) must be modified accordingly. 

The rigid body modes associated with shape function must be eliminated in order to define elastic 
deformations uniquely. This can be achieved using a specified body frame, a set of mode shapes, or an 
imposed set of boundary conditions. These three approaches are essentially the same. For example, im- 
posing the condition that the body frame of a beam be rigidly attached to a node at one end of the beam 
is equivalent to selecting a set of cantilever modes. In this formulation, it is assumed that a body frame 
is rigidly attached to a node associated with the body. This is not a limitation of this formulation. If an- 
other set of body frames is more appropriate or if another set of modes are more desirable, then one can 
develop a proper transformation relation between nodal elastic variables and minimum elastic coordinates 
of the body. The transformation matrices appearing in such relations are in general constant and linear, 
and therefore such transformations can be applied throughout the formulation without any change in the 
formulation presented in this report. 


3. VIRTUAL WORK DUE TO ELASTIC FORCES 
The virtual work due to elastic forces in body i is given by (refs. 7, 8) 

8Wl= f. Se‘VW (3.1) 

Jv i 

where e‘ and a* are, respectively, vectors of strain and stress at a point in body i. The virtual work due to 
elastic forces considered in equation (3.1) is general and includes nonlinear stress-strain relationship and 
the viscoelastic materials. For simplicity, only a linear stress-strain relation is considered hereafter. Using 
Hooke’s law, the stress-strain relation can be written as 
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<7 l = Et* 


(3.2) 


and for small elastic deformations, the strain displacement relation can be written as 


e’ = Du (3.3) 

where E is the modulus of elasticity matrix and D is a differential operator. Using equation (2.3), the strain 
displacement relation can be expressed as 


e’ = (3.4) 

Substituting equations (3.2) and (3.4) into equation (3.1) yields 

8Wi = ( 6q {r (DN i ) T E(DN i )q i dv = 8q <r K[ e q { (3.5) 

Jv 1 

where K * e is the elastic stiffness matrix of body i. The stiffness matrix K\ e is same as the finite-element 
stiffness matrix appearing in standard linear structural analysis and is a constant matrix. Therefore, one 
needs to generate this matrix only once. Using the generalized coordinate vector P\ the virtual work 8W' e 
is 


8W' e = 8P' T K i P i 


where 


K' = 


0 0 0 
0 0 0 
0 0 Ki e 


(3.6) 


(3.7) 


There is no stiffness coupling between rigid body and elastic coordinates because the strain depends entirely 
on elastic deformations. 

The total virtual work due to elastic forces is given by 


8W e = Y, 8w< e = Z) 8 >P <r K i P i 

i i 


(3.8) 
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4. VIRTUAL WORK PRODUCED BY EXTERNAL FORCES 


The external forces considered here may be point loads, body forces, or surface forces, or a combina- 
tion of all three. Body forces include the portion of the load attributed to gravity, surface forces may be 
aerodynamic drag and lift and other surface tractions, and point loads may arise when certain control ele- 
ments are present. Furthermore, forces due to several other elements, such as translational and/or rotational 
spring-damper actuators, can also be considered. 

Let F£ and F, be, respectively, the body force and the surface force acting on body i. The virtual work 
produced by these forces is given by 


8Wf = f Ff ■ 8rdv' + [ F? • 8rds { (4.1) 

Jv Js 

In the case of point loads, force F£ is expressed using a delta function. Substituting the value of the 
virtual displacement from equation (2.13) into equation (4.1), one finds 

8Wi = |jf Ff [ I B\6, s) A'N' ] dv { + J* F? [ I B { (8, s ) A'N' ] ds‘‘] • 8P i 

= 8P {r -Q i (4.2) 

where 



* I * 


' I 

b <t 

Ft-dv*+ [ 

B? 

. N ,r A‘ r . 

Js 

_ N' T A< r . 


Fi ■ ds< 


(4.3) 


is the vector of generalized forces. In equation (4.3), the first row represents the vector of translational 
force along the inertial axis, the second row is associated with the torque about the body axis, and the third 
row is a vector of generalized forces associated with virtual elastic deformations. The torque about the 
body axes, which is constant in rigid body analysis, is no longer constant due to elastic deformations. The 
total virtual work due to all external forces acting on the system is given by 


8w t = 5W t = Z) spiT • <2‘ 

t t 


(4.4) 


The virtual work due to the workless constraint forces does not appear in the above expression. 
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5. FORMULATION METHODS AND KINEMATIC CONSIDERATIONS 


Having obtained the virtual work expressions due to inertia forces, internal elastic forces, and external 
forces, one can proceed to derive the equations of motion of the system. For dynamic equilibrium of the 
system, the sum of the virtual work due to inertia and elastic forces must be equal to the virtual work due 
to external forces. This implies 


8W I + 6W e = SWi (5.1) 

Substituting the values of 8Wi, 8W e , and 8Wi from equations (2.28), (3.8), and (4.4), leads to 


NB 

E & 

i=l 


[M i P i + F' + K i P i - Q { 


= 0 


(5.2) 


If the virtual displacements 8P' ( i = 1, ..., NB) are independent of each other, then equation (5.2) is 
satisfied only if the coefficients of 8P ‘ are zero. This implies 


M i P i + F i + K i P i - Q* = 0 (i=l,...,NB) (5.3) 

Equation (5.3), however, is not true if the components are kinematically constrained. For example, in 
a helicopter, the fuselage and a rotor are kinematically constrained to allow only one relative degree of 
freedom with respect to each other. Similarly, one may have a universal joint between a rotor blade and 
the rotor hub. Such kinematic joints impose nonlinear algebraic constraints that must be satisfied. This 
implies that the generalized coordinates P’(i = 1 , ..., NB) are no longer independent, and the coefficient 
of 8P ‘ in equation (5.2) can not be set equal to zero. Therefore, equation (5.2) must be modified to account 
for the kinematic constraints. 

The kinematic constraints between adjacent bodies can be incorporated into the equations of motion 
using one of the following three approaches: 1) the Lagrange multiplier approach, 2) explicit reduction of 
the equations of motion in terms of independent generalized coordinates, and 3) the hybrid approach. 

In the first approach, a set of algebraic constraints that defines the kinematics of the system is devel- 
oped. This set is then introduced into equation (5.2) using the Lagrange multiplier approach. This leads 
to a system of differential and algebraic equations that must be solved simultaneously. The structure of 
the algebraic constraint equations depends on the type of the kinematic constraints (or joints) and is to- 
tally independent of the type of system considered. Furthermore, the algebraic equations associated with 
a kinematic joint are functions of the generalized coordinates of only the two components that are partici- 
pating in the kinematic joint. Therefore, in this approach, a library of constraints can be developed and the 
constraints between two adjacent elements can be imposed in an automated fashion. 
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In the second approach, the differential equations of motion are developed in terms of a set of inde- 
pendent coordinates. In most cases it is difficult to express all virtual work expressions entirely in terms 
of independent coordinates. Therefore, the dependent virtual displacements and the dependent acceler- 
ation in equation (5.2) are replaced, respectively, by independent virtual displacements and independent 
accelerations, and the mass matrix and other matrices are allowed to be functions of dependent and inde- 
pendent coordinates. This is always possible because virtual displacements are linearly related. Similarly, 
the components of the acceleration vector are also linearly related. One advantage of this approach is that 
the dynamics of the system are expressed in terms of the minimum number of second-order differential 
equations. In this approach, however, the sparsity of the matrices is lost and several submatrices which 
were constant at the component level do not remain constant in the assembled system of equations. Fur- 
thermore, this approach yields special purpose codes of hardwired mathematical models that are generally 
difficult to update. 

In the third approach, the system is divided into a number of groups consisting of one or more com- 
ponents in each group such that no component is a member of two or more groups. The virtual work 
expressions are then developed for each group in terms of its minimum number of generalized coordi- 
nates. Then each group is dynamically coupled to the adjacent group using Lagrange multipliers. Such 
groups are called superelements. This approach attempts to take advantage of the above two approaches. 
A set of superelements that accounts for every kinematic joint is not yet available. Furthermore, such code 
can be used efficiently by the most experienced users only. In this report, only the first approach will be 
expanded. 


6. CONSTRAINT EQUATIONS 


In this section, a set of constraint equations for a revolute joint that includes the coupling between a 
fuselage and a rotor is developed. This is done to illustrate how a particular set of constraints is derived and 
to show by implication that a library of constraint cases can be derived and treated as constraint elements 
in the modeling processes. 

Figure 2 shows a simplified representation of a fuselage (body i) and a rotor (body j). Let the two 
points P and Q lie along the axis of the rotor, and be common to the two bodies. The kinematic compati- 
bility condition between the two bodies requires that: (1) the point P on body i and on body j should be 
connected, and (2) the vector PQ in body i and in body j should always be collinear. Mathematically, the 
first condition implies 


R + A\ Sopi + N\ P) • < 7 ‘) - R - A’(s opj + N’(P)q’) = 0 (6.1) 

where Sopi and Sopj are, respectively, the undeformed position vector of point P with respect to the global 
axes of body i and body j, and N'(P) and N>(P) are, respectively, the value of the shape function N 1 
and at point P. 
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The vector PQ in body i is given by 


( PQ) *' = A\ ( Soqi - Sopi ) + ( N\ Q) - NXP ) ) q*) (6.2) 


Similarly, the vector (PQ) in body; is 


( PQY = - Sop; - ) + (AP'(Q) - NXP))q j ) (6.3) 


The second compatibility condition then requires that 


(PQ)' x (PQ) ; = 0 (6.4) 

Equation (6.4) gives three constraint conditions; however, only two of them are independent. Any two 
conditions from equation (6.4) are sufficient to impose the second compatibility condition. The five con- 
straint equations, (three from equation (6.1) and two from equation (6.4)) define the kinematic constraint 
between a fuselage and a rotor. This approach allows the elastic deformations at the kinematic joint. If 
body i and/or body j are rigid then the corresponding constraint equations are obtained by setting the elastic 
deformation of the rigid body equal to zero. Also points P and Q need not be node points. However, if 
points P and Q are chosen at two nodes, then equations (6.1) and (6.4) are further simplified. 

The constraint equations developed above are general and are applicable to all revolute joints ( and 
therefore it applies to a fuselage and a rotor pair also). Furthermore, these constraints are functions of the 
generalized coordinates of only the two bodies. Therefore, only one set of constraint equations for each 
revolute joint is required. 

The set of constraint equations developed above is known as a revolute joint constraint element. A 
helicopter system may include other types of kinematic joints, such as universal joints (e.g., a universal 
joint between a rotor and a blade), spherical joints, translational joints, cylindrical joints, and spherical- 
cylindrical joints. Consideration of these joints and several other kinematic joints can be quite important 
in modeling a helicopter. Constraint elements for these joints can be developed from purely geometri- 
cal considerations. For example, a constraint element for a universal joint, as shown in figure 3, can be 
developed by imposing these conditions (1) point P in body i and in body j should be connected, and 
(2) the vector PA in body i should be normal to the vector PB in body j. Constraint equations developed 
in this manner can be incorporated in a library of nonlinear constraint elements. These elements then can 
be used to set up automated kinematic constraints between two adjacent bodies. 

The complete set of constraint equations will be represented by 


<KP,t) = o 


(6.5) 
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where $ is a m x 1 vector of constraints, P is a composite vector of P* given by 

P=[P' T ... P NBT ] T (6.6) 

and t is the time parameter. The time parameter has been included here to emphasize that the constraints 
may explicitly depend on time. 


7. SYSTEM EQUATIONS OF MOTION 


In deriving differential equations of motion via the virtual work approach (eq. (5.2)), the requirement 
for constraints (eq. (6.5)) does not appear until the last step, when the variations of 8P i are considered 
independent of each other. Recall that the virtual displacements 8P ' in equation (5.2) are not independent 
of each other because of the presence of constraints. Therefore the coefficients of 6P* in equation (5.2) 
cannot-esettozero. Inordertodevelopthevirtualworkexpressionintermsofindependentvirtual 
displacements, it was noted that equation (6.5) is equivalent to a differential equation, 


X.w ip,+d i^° 


(7.1) 


By definition, the virtual displacement is one in which the time is considered frozen, hence the virtual 
displacements 8P * must satisfy 




If equation (7.2) is true, then it is also true that 


(7.2) 



(7.3) 


where A is a to x 1 vector of unknown Lagrange multipliers. Equation (7.3) can be adjoined with equation 
(5.2), yielding the relation 


NB 


M ipi + pi + pipi _ Qi + ( 



= 0 


(7.4) 


Note that 6P‘(t = 1 , ..., NB) are still not independent. Following an approach presented in reference 9, 
it can be shown that 
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(7.5) 


M ipi + F i + K ipi _ Q< + ( ^) T \ = 0 

Equation (7.5) is the desired set of differential equations of motion. If the idealized system has n number 
of generalized coordinates and m Lagrange multipliers, then equation (7.5) together with equation (6.5) 
provides (n+ m) equations in (n+ m) unknowns. Equations (7.5) and (6.5) together define a set of 
differential and algebraic equations that must be solved simultaneously. 

In the present formulation, the matrices M\ F\ K\ and Q* remain completely unaffected by the 
kinematic joints appearing in the system. As a result, the constant submatrices appearing in matrices M\ 
F\ K\ and Q‘ remain constant and the sparsity of these matrices is preserved throughout the dynamic 
simulation. Further, the components of the system can be joined to adjacent elements by selecting proper 
constraints from the library of constraint elements. This formulation also allows generation of the equations 
of motion of the system with any number of flexible and/or rigid components. An additional benefit is that 
the vector 


( 



(7.6) 


provides reaction forces at the kinematic joints. 

The derivation of equations of motion presented is general, and one can consider any set of rotational 
degrees of freedom and arbitrary finite-element shape functions. The specific form of the various matrices 
appearing in the above formulation will depend on the rotational degrees of freedom and the finite-element 
shape functions considered. 


8. NUMERICAL COMPUTATION METHODS 


The mathematical formulation presented in the previous sections leads to a set of differential and al- 
gebraic equations. In this section, a brief review of numerical methods used to solve the set of differential 
and algebraic equations is presented (eqs. 6.5 and 7.5). The numerical schemes considered here are those 
which have been used extensively in multibody dynamics. These numerical schemes can be divided into 
the following three groups: schemes which satisfy the constraints explicitly at every step; schemes which 
attempt to satisfy the constraints implicitly; and hybrid schemes, in which the above two approaches are 
combined. 

Numerical schemes that satisfy the constraints explicitly at every step are as follows: lower upper (LU) 
factorization technique (LUFT) (ref. 10), singular value decomposition technique (SVDT) (refs. 11, 12), 
QR decomposition technique (QRDT) (ref. 13), null space method (NSM) (refs. 14, 15), and tangent 
coordinate method (TCM) (ref. 16). 

In each of these schemes, the independent coordinates are integrated using a differential-equation so- 
lution subroutine, and dependent coordinates are obtained using an iteration scheme such as the Newton- 
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Raphson scheme. In LUFT, the independent coordinates are identified using LU factorizations of the 
constraint Jacobian matrix. The theory of differential equations in manifold shows that the solution of the 
dynamical equations lies on the constraint surface. In SVDT, QRDT, and TCM the tangent coordinates are 
generated using, respectively, Singular Value Decomposition, QR Decomposition, and Gramm-Schmidt 
Orthonormalization of the Jacobian matrix. 

Schemes that attempt to satisfy the constraints implicitly are: the constraint stabilization technique 
(CST) (refs. 17, 18), and the penalty function method (PFM) (ref. 19). Some variations of these schemes 
have also been developed (ref. 20). In these techniques, attempts are made to implicitly satisfy the con- 
straints. This is achieved by introducing artificial damping into the system differential equations of motion. 

In the hybrid schemes, the Constraint Stabilization Technique is used to advance the integration from 
one step to the next step and the constraints are checked at every step. If the constraint violations goes above 
the error tolerance, then the dependent coordinates are corrected to satisfy the constraints and integration 
is continued. 

It should be mentioned here that the integration schemes considered are those that can be completely 
isolated from the mathematical formulation of the equations of motion. This is quite important for further 
extension of the code. 


9. CONCLUDING REMARKS 


A general mathematical formulation for automated computer-aided analysis of a helicopter system, 
whose components may undergo large translations and large rotations, has been presented. In this formu- 
lation, a helicopter system is considered as a collection of several components (bodies) interconnected by 
kinematic joints to allow the kinematic degrees of freedom between two adjacent components. A body 
reference frame is associated with each component. Two sets of generalized coordinates, namely rigid 
body and elastic, are used to locate a point on the component. The rigid body coordinates define the po- 
sition and the orientation of the body frame. Each flexible component is idealized as a collection of finite 
elements rigidly connected at the nodal points. The elastic deformations at the nodal points in a body are 
defined with respect to the body frame and a transformation matrix from the body frame to the inertial 
frame is used to define the elastic deformations in the inertial frame. The deformations between two nodes 
are approximated using shape functions that are defined in the body frame. The generalized coordinates 
defined above are used to develop expressions for virtual work produced by inertia forces, elastic forces, 
and external forces. The kinematic constraint conditions between two adjacent bodies are expressed by a 
set of nonlinear constraint equations. The virtual work approach is used to derive the equations of motion 
and the Lagrange multiplier approach is used to adjoin the constraint equations to the equations of motion. 
The model thus combines the traditional finite element approaches that are capable of handling small elas- 
tic deformations with the multibody approaches that allow its components to undergo large translations 
and large rotations. 

In the mathematical development, the role of a body axis needs to be emphasized. In this formulation, 
the elastic deformations between the nodes are approximated using shape functions. This is possible pri- 
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manly because the strain and the rotation of an element about its axis is small (because of the choice of 
body frame), and the elastic coordinates are defined with respect to the body frame. 

Three features of the transformation matrices that map the body coordinates to inertial coordinates 
are noteworthy. First, the transformation matrices are used to define the configuration of each component 
in a single inertial frame. This leads to a single basis-frame in contrast to a vector approach where dif- 
ferent vectors may have different basis-frames. Second, the kinematic relations are developed by direct 
differentiation of the configuration field without consideration of the coordinate systems used to define the 
parameters. This is significantly different from a vector approach where the total derivative and the deriva- 
tive with respect to a rotating body axis are not the same. Third, by using transformation matrices, a set of 
rotational coordinates can be selected that is completely isolated from the other generalized coordinates. 
The transformation matrices therefore provide flexibility in the choice of rotational coordinates without 
altering the formulation. 

In this model, the generalized coordinates of a component are defined independently of the generalized 
coordinates of the other components. As a result, a generic component is used to develop a set of mathe- 
matical models for inertia effects, internal elastic effects, and effects due to external forces. Each model 
can be developed in separate modules that can be used for each component in the system. This model of 
a component has two further advantages. First, the constant matrices associated with a component remain 
constant, and second, the sparsity of the component matrices is preserved. 

In this formulation, two components are connected in an arbitrary fashion. This connection is accom- 
plished eveloping a set of nonlinear constraint equations and using a Lagrange multiplier technique. The 
constraint equations are generated purely from geometrical considerations. A methodology of generating 
nonlinear constraints for revolute joints and universal joints have been presented. Nonlinear kinematic 
constraints associated with other types of kinematic couplings, such as articulated joints, spherical joints, 
translational joints, and other joints can be developed and incorporated as a library of constraint elements. 
This library then can be used to enforce the kinematic constraints in an automated fashion. The Lagrange 
multipliers yield the constraint forces at the joints. 

Finally, the two most important features of this formulations are the generality of the mathematical 
model and the modularity of the code. Since this model is general, it allows one to model a helicopter 
system with arbitrary topological configurations. The generality of this approach allows for theoretical 
growth of this model without complete reformulation of the problem. This capability is very important 
when other significant physical aspects of helicopter dynamics must be incorporated in the software as the 
technology grows and the demand for high accuracy becomes important. The modularity of the code in 
this approach lends itself to further enhancements. As a result, each module in the code can be developed, 
checked, and modified independently. 
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